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Abstract 

A detailed study is carried out for the relativistic theory of viscoelasticity which 
was recently constructed on the basis of Onsager's linear nonequilibrium thermody- 
namics. After rederiving the theory using a local argument with the entropy current, 
we show that this theory universally reduces to the standard relativistic Navier-Stokes 
fluid mechanics in the long time limit. Since effects of elasticity are taken into account, 
the dynamics at short time scales is modified from that given by the Navier-Stokes 
equations, so that acausal problems intrinsic to relativistic Navier-Stokes fluids are sig- 
nificantly remedied. We in particular show that the wave equations for the propagation 
of disturbance around a hydrostatic equilibrium in Minkowski spacetime become sym- 
metric hyperbolic for some range of parameters, so that the model is free of acausality 
problems. This observation suggests that the relativistic viscoelastic model with such 
parameters can be regarded as a causal completion of relativistic Navier-Stokes fluid 
mechanics. By adjusting parameters to various values, this theory can treat a wide 
variety of materials including elastic materials, Maxwell materials, Kelvin- Voigt ma- 
terials, and (a nonlinearly generalized version of) simplified Israel-Stewart fluids, and 
thus we expect the theory to be the most universal description of single-component 
relativistic continuum materials. We also show that the presence of strains and the 
corresponding change in temperature are naturally unified through the Tolman law in 
a generally covariant description of continuum mechanics. 
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1. Introduction 



The dynamics of fluids at large scales is universally described by the Navier-Stokes equations, 
which represent the regression to a global equilibrium with transfers of conserved quantities 
(such as energy- momentum and particle number) among fluid particles PQ. This can be 
formulated in a generally covariant way, but it is known that there arises a problem of 
acausality. In fact, the obtained equations for the propagation of disturbance are basically 
parabolic and thus predict infinitely large speed of propagation for infinitely high frequency 
modes, leaving light cones. One should note here that this does not imply the breakdown of 
the internal consistency of the description because the Navier-Stokes equations are simply an 
effective description at large spacetime scales and need not describe high frequency modes 
correctly. However, this is still troublesome when adopting the equations in numerical 
simulations; the initial value problems are ill posed, and unacceptable numerical solutions 
can be obtained easily. 

To remedy the problem, Miiller, Israel, and Stewart [21 El S] extended the theory by 
treating the dissipative part of stress tensor, , and the heat flux (for the Eckart 
frame) or the particle diffusion current (for the Landau-Lifshitz frame) as additional 
thermodynamic variables on which the entropy density can depend. This prescription is 
based on the so-called extended thermodynamics and corresponds to taking into account 
higher derivative corrections to the effective theory. It has been shown that such modified 
theories have a good causal behavior and that linear perturbations around a hydrostatic 
equilibrium obey hyperbolic differential equations. This is now regarded as a fundamental 
framework for the numerical study of relativistic viscous fluids. 

Meanwhile, modifications of the Navier-Stokes equations have also been studied in the 
area of rheology, and the materials treated there are generically called viscoelastic materials 
or viscoelastic fluids. Historically, viscoelasticity was defined by Maxwell in the 19th century 
as the characteristic property of such continuum materials that behave as elastic solids at 
short time scales and as viscous fluids at long time scales [HE]- In 1948, Eckart proposed 
a theory of elasticity and anelasticity [6], which describes the nonrelativistic dynamics of 
single-component viscoelastic materials and was reinvented recently [7J in the light of the 
covariance under foliation preserving diffeomorphisms. In this description, elastic strains 
(or equivalently, the "intrinsic metric" defined below) are introduced as additional thermo- 
dynamic variables, as in the theory of elasticity. As explicitly shown in [7J, this theory of 
viscoelasticity contains the theory of elasticity and the theory of fluids as special limiting 
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cases, and correctly reproduces the Navier-Stokes equations in the fluid limit. Furthermore, 
as was pointed out in [5], since the dynamics at short time scales is dominated by elastic- 
ity, shear modes of linear perturbations around a hydrostatic equilibrium obey differential 
equations with second-order time derivatives (in contrast to the equations obtained from 
the Navier-Stokes equations that contain only a first-order time derivative), so that causal 
behaviors for large frequencies are significantly improved. 

Recently, on the basis of Onsager's linear regression theory on nonequilibrium thermo- 
dynamics [HI dni EH C2] , the present authors proposed a relativistic theory of viscoelasticity 
[T3] which generalizes the theory of elasticity and anelasticity [HI E] in a generally covariant 
form. In the present paper, after rederiving the theory relying on a local argument with 
the entropy current, we study the detailed properties of relativistic viscoelasticity. We show 
that fluidity is universally realized in the long time limit and also that acausal problems 
disappear for a wide region of parameters. Thus, the relativistic theory of viscoelasticity 
with such parameters can be regarded as a causal completion of relativistic Navier-Stokes 
fluid mechanics, and we expect that it could be used as another basis in the numerical study 
of relativistic viscous fluids. 

This paper is organized as follows. In Sec. |2] we rederive the viscoelastic model of [13] 
using a local argument with the entropy current. We also show that the presence of strains 
and the corresponding change in temperature are naturally unified in a generally covariant 
description of continuum mechanics. In Sec. [3] we consider the long and short time limits of 
our viscoelastic model. We prove that the model universally gives relativistic Navier-Stokes 
fluids in the long time limit. In Sec. H]we show that when some parameters take specific 
values, our viscoelastic model reduces to (a higher-dimensional extension of) the nonlin- 
ear generalization of the simplified Israel-Stewart model [2]. In Sec. E]we consider linear 
perturbations around a hydrostatic equilibrium in Minkowski spacetime. The dispersion 
relations show that the evolutions are certainly stable. Although the wave equations for the 
linear perturbations are not always hyperbolic, if some parameters are chosen appropriately 
(including the parametrizations for the simplified Israel- Stewart model) they become sym- 
metric hyperbolic and thus free of acausality problems. Section [6] is devoted to conclusion 
and discussions. 
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2. Relativist ic viscoelastic mechanics 



In this section, we rederive the fundamental equations for relativistic viscoelastic mechanics 
using a local argument with the entropy current. In Appendix |A] we show that the present 
formulation is equivalent to the "entropic formulation" proposed in our previous paper [TH] 
which is based on Onsager's linear regression theory. 

2.1. Definitions 

We start by giving a brief review on the generally covariant definitions of viscoelastic mate- 
rials [T3] . 

Geometrical setup 

We consider a single-component continuum material living in a (D + l)-dimensional 
Lorentzian manifold M. . The local coordinates are denoted by x M (// = 0, 1, • • • , D) , and 
the background Lorentzian metric with signature (— , +,•••,+) by g^x) . Following the 
convention of Landau and Lifshitz [Tj, we define the velocity field u = u fl (x)d tl from the 
momentum (D + l)-vector p^ as 

u"(x) = g^(x) p„(x)/e(x) = p' I (x)/e(x) , (2.1) 

where e(x) = \/—g^ u (x) p p (x) p u (x) is the proper energy density. Note that u^(x) is normal- 
ized as g fJiU (x)u' J '(x)u u (x) = — 1 . Here and hereafter indices are subscripted (or superscripted) 
always with g py (or with its inverse g^). 

Assuming that the velocity field is hypersurface orthogonal, we introduce a foliation of 
Ai consisting of spatial hypersurfaces (timeslices) orthogonal to u M . We parametrize the 
timeslices with a real parameter t and denote them by E t . We exclusively (except for Sec. 
[5]) use a coordinate system x = (x^ 1 ) = (x°, x) such that x° = t , x = (x l ) (i = 1 , • • • , D) , 
for which the shape of the material at time t is given by the induced metric on : 

h^{x) = g^u(x) + u^x) u v {x) = g^(x) + M^^i^j . (2.2) 

e [x) 

We also define the extrinsic curvature of the hypersurface as half the Lie derivative of 
h pu with respect to the velocity field u = u^d^ : 

K„ v = - £ u V = \ hfhj (V p u a + V a u p ) . (2.3) 
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This measures the rate of change in the induced metric as material particles flow along 
u M . Note that this tensor is symmetric and orthogonal to u^, K^ u u v = . 

In the ADM parametrization, the metric and the velocity are represented with the lapse 
N(x) and the shifts N*(x) (i = 1, ■■• ,D) as 

ds 2 = g^{x) dx^ dx v = -N 2 (x) dt 2 + h^x) (dx i - N*(x) dt) (dx j - N j (x) dt) , (2.4) 
1 NUr) 

u = u»(x) d, = — - d + — di u,(x) dx» = -N(x) dt) . (2.5) 

The volume element on the hypersurface is given by the D-form \/hd D x = y/det(hij) d D x = 
N~ l ^jd D x. 

With a given foliation, we still have the symmetry of foliation preserving diffeomorphisms 
that give rise to transformations only among the points on each timeslice. Using this residual 
gauge symmetry we can impose the synchronized gauge, N l (x) = , so that the background 
metric and the velocity field are expressed as 



u = u^x) 8^ = T77 - K — = — , (2.7) 



ds 2 = g^(x) dx 11 dx u = -N 2 (x) dt 2 + h^x) dx i dx j , (2.6) 
1 d d_ 
N{x)di ~ ~{fr 

where r is the local proper time defined by dr = N dt . In this gauge, due to the relation 
d/dt = N(x)d/dr, the proper energy density e(x) measured with the proper time r is 
related to the energy density e(x) measured with time t as 

e(x) = N(x) e(x) . (2.8) 

Note that e(x) includes the gravitational potential through the factor N(x) . Accordingly, 
the local temperature T measured with r is related to the temperature T measured with t 
through the following Tolman law: 

T(x) = N(x)T(x) . (2.9) 



Definition of (relativistic) viscoelastic materials 

According to the definition of Maxwell, viscoelastic materials behave as elastic solids at 
short time scales and as viscous fluids at long time scales (see, e.g., Sec. 36 in [5]). In order 
to understand how such materials evolve in time, we consider a material consisting of many 
molecules bonding each other and assume that the molecules first stay at their equilibrium 
positions in the absence of strains (as in the leftmost illustration of Fig. [Tj) [3 [8] . We now 
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deformation 
>- 




Figure 1: Processes of deformation and stress relaxation 



suppose that an external force is applied to deform the material. An internal strain is then 
produced in the body, and according to the definition, the accompanied internal stress can 
be treated as an elastic force at least during short intervals of time. However, if we keep the 
deformation much longer than the relaxation times (characteristic to each material), then 
the bonding structure changes to maximize the entropy, and the internal strain vanishes 
eventually as in the rightmost of Fig. [TJ The point is that two figures (the central and the 
rightmost) have the same shape (same induced metric) h^ u , but different bonding structures. 

The internal bonding structure can be specified by the intrinsic metric , which mea- 
sures the shape that the material would take when all the internal strains are removed 
virtually [61 [TJ . For the example given in Fig. (TJ the intrinsic metric for the center illustra- 
tion is given by the induced metric for the leftmost illustration, while the intrinsic metric for 
the rightmost illustration agrees with the induced metric for itself. Thus, the plastic (i.e., 
nonelastic) deformation from the center illustration to the rightmost illustration is described 
as the evolution of the intrinsic metric^] 

Its generally covariant generalization can be defined in the following way. Suppose that 
we have two adjacent, spatially separated spacetime points P and Q , each of which rep- 
resents a point on the trajectory of a material particle (see Fig. [2]). By denoting their 
coordinates by x = (x^) and x + dx = (x M + dx^) , respectively, the distance between P and 
Q in the real configuration is of course given with the metric as (the square root of) 



We now virtually remove all the strains in a sufficiently small spacetime region including the 
two points. Then P and Q would move to other positions P and Q, whose coordinates we 
denote by x = (a; M ) and x + dx = (x M + dx M ) , respectively. This correspondence defines a 
local map / : x i— >• x — x(x) , with which we define the intrinsic metric g^x) as the metric 

h^v is also called the "strain metric" and was first introduced by Eckart to embody "the principle of 
relaxability-in-thc-small" in anelasticity [5] . Some examples of the explicit form of h^ v and h^ v under various 
deformations can be found in [71 [8] . 



ds 2 = g^ u (x) dx M dx 



(2.10) 
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x = x{t) 




Figure 2: Real (x M (r)) and virtual (x M (r)) trajectories of material particles. The 
distance between P and Q gives the definition of the intrinsic metric . 



measuring the virtual distance between P and Q (or, as the pullback of the metric g^ u for 
the map; g^ v = /*flf^)fi 

ds 2 = g pa (x) dx p dx° = g per (x{x)) — — dx» dx v 

= g^ix) dx^ dx u . 

With the velocity vector u = u^(x) , we parametrize g^ u as 



(2-11] 



9i* 



1 + 26) Ufj,u v - e p u u - e u u M + h 

V 



[1 I 10 flV 



(e/ = 0, h^u u = 0, h^u^O). 



(2.12) 



The strain tensor is then introduced as 



1 



where 



Epy{x) = - (g^(x) - g^(x)) 

= 9 u p u u + - [e^Uy + e v u M ) + 



■<iw(b) = - (h^(x) - h pv (x)) 



1 



(2.13) 



(2.14) 



is the spatial strain tensor. Note that if we define the extrinsic curvature associated with 
the spatial intrinsic metric h pu as 



K n» = o £ u h ^u = ^ {u x d x h pu + dpU X h Xu + d u u x h pX ) 



the following identity holds: 



£ e = K - K 



(2.15) 



(2.16) 



2 As in the standard theory of elasticity [5], there may be an arbitrariness in defining x M , but the intrinsic 
metric can still be defined uniquely. 
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A viscoelastic material is a thermodynamic system consisting of material particles as its 
subsystems. While the system regresses to a thermodynamic equilibrium, one can imag- 
ine that the virtual trajectory of each material particle approaches its real trajectory, so 
that the strain tensor approaches zero. Such an irreversible process is called plas- 
tic (i.e., nonelastic), and thus we see that the dynamics of includes plastic evolutions 
(in addition to reversible, elastic evolutions). In the following discussions, we assume that 
E^ v = [e^y ,6^,9) are all small quantities, such that their nonlinear effects can be neglected. 
We shall denote the contraction of a spatial tensoio with g^ v by tr A , so that 

tre = g^e^ = h^e^ , tvK = g^K^ = h^K^ . (2.17) 

We close this section by explaining the physical meaning of the strain tensor E^ u = 
{tfiv ? £(i ? 0) ■ The spatial strain tensor stands for the standard strains, measuring the 
difference between the induced metric h^ u and the spatial induced metric . One can 
easily see that the quantity represents the relative velocity of a material particle in its 
real trajectory with respect to that in its virtual trajectory, = u^—u^ 1 = dx^/dr — dx^/dr , 
where r is a common proper time (see Fig. [2]). In order to understand the meaning of 9, 
we first recall that the covariant vector is expressed as dx M = —N dx° . We can then 
rewrite ds 2 and ds 2 as 

ds 2 = - N 2 {x) (dx ) 2 + h^(x) dx" dx u , (2.18) 
ds 2 = - (1 + 29{x)) N 2 (x) (dx ) 2 + 2N(x) e^x) dx" dx° + h^(x) dx" dx v , (2.19) 

with h^ v dx" dx u = hij (dx 1 — N l dx°) (dx- 7 — A" dx°) and a similar (but a bit more com- 
plicated) expression for 2iV dx" dx° + dx" dx v . These equations mean that N = 
yl + 29 N ~ (1 + 9) N represents the lapse function for the intrinsic metric. Then, through 
the Tolman law, we can relate the virtual temperature T observed in the absence of strains 
to the actual temperature T as JVT = NT(= T) . We thus obtain the relation 9 = 
(N 2 /N 2 - l)/2 = (T 2 /T 2 - l)/2 ~ (T - f )/f , and conclude that the scalar 9 expresses 
the increase of the temperature due to strains. This conclusion shows that the presence 
of strains and the corresponding change in temperature are naturally unified in a generally 
covariant description of continuum mechanics. 



3 By spatial we mean that is orthogonal to u M , A^ u u 1 ' = = A^u^ . Recall that = —u^u v -\-h ,lv 
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2.2. Entropy production rate 



As was adopted in [13], in order to develop thermodynamics in a generally covariant manner, 
it is convenient to distinguish density quantities from other intensive quantities, and, by 
multiplying them with the spatial volume element y/h , we construct new quantities which 
are spatial densities on each timeslice. For example, the entropy density s, the energy- 
momentum density p M , and the number density n are density quantities, and for them we 
introduce the following spatial densities: 

s = y/hs, p^ = y/hp^, n = yhn. (2.20) 

We assume that each material particle is in its local thermodynamic equilibrium, and that 
the local entropy s is a function of p^ , n , and g^ as well as of the strain tensor Ep V = 

{^flU -i £(M -i 9) ■ 

s(x) = s(E^ u (x),p^(x),h(x), g^ u (x)) . (2.21) 

We further assume that s depends on p^ only through the local proper energy e(p M , g^ v ) = 
yf—g^ PuPu , so that s can also be expressed as 

s(x) = a(E^ u (x) , e(x) , n(x) , fiv( x )) 

= v(e,u,(x) , e M (sc) , 9{x) , e(p^(x), g^x)) , h(x) , g^(x)) . (2.22) 

Since we are only interested in linear nonequilibrium thermodynamics, we only need to 
expand s in E^ to second orderQ 

s = (terms independent of E^ v ) 

- — [2X^)8^ + A 2 + 7i (tr ef + 2 72 (tr e) 9 + 73 9 2 ] . (2.23) 

We require the stability of the system under the change in strains E^ v , so that the constants 
Ai and A 2 are non-negative, and the matrix 7 = (72 73 ) is positive semidefmite. 

Then the fundamental thermodynamic relation can be written as 

S3 = ~ Y S P» ~ f ^ + ^ T (q) 89^V 

~2\ x e<""> 8e { , v) - ^ ( 7l tr e + l2 9) 5(tr e) 



'For a tensor , we define A (fll/) = (1/2) h£ hf [A pa + A ap - (2/D) h a ? A a/3 h pa ] 
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(2.24) 



rj-\ ^ ft rj~i 

Here the temperature T , the chemical potential /i and the quasiconservative part of the 



stress tensor, , are defined as 



da_ 

de 



1 



9a 
<9n 



<9cr v/i 



T ' %,"2T r S' (2 ' 25) 

where we require that be orthogonal to m m , u v = . The quasiconservative part of 
the energy-momentum tensor is then defined as 



In deriving Eq. ( 12.24(1 . we have used the relations 



-u 



— _ I* it 



We now set the variation in Eq. (" 12.24(1 to be 5 = £ u . We then obtain 



(2.26) 



(2.27) 



Vh 



- y V m (r, U ") - f V> «") + - 7^ ^ 

2A 1 

- ~y £<""> £ U £ (H - - (71 tr £ + 72 0) £„(tr e) 

- -^T £ M - ^ (73 + 72 tr e) £ u 9 



(2.28) 



Here we have used the identities for Lie derivatives: 

£ u s = ^V M (><), £ u p v = VhV li (p„u ti ), £ u n = VhV fl (nu^) J (2.29) 

which can be shown by using the identities £ u \fh = y/hV ^u 11 and p^V^u^ = 0. Note 
that tr(£ M e M1 ,) = h^ v £ v e iW can be replaced by £ u (tre) in our approximation because the 
difference £ u (tre) — ti[£ u e^ v ) = (£ u h fJ,v ) = —2K^ u e^ u is of higher orders. 

The full energy-momentum tensor T^ v and the full number current n M are given by 

= e u»u v + , n^nu^/, 

= = (2.30) 

where t^ u and are the stress tensor and the diffusion current, respectively. Then, by 
introducing the entropy current 
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(2.31) 



5 We here use a convention that the quasiconservative stress tensor does not include stresses originated 



from strains. 



and by using Eq. ( I2.28P together with the current conservation laws 



V V T» V = , vx = 



the local entropy production rate can be evaluated as 

1 



(2.32) 



£ u t 



T v (q) 
2Ai 

~Y 

)l 
r 

(ju>) 



)K, v +^d,{-^ 



- — (71 tr e + 72 9) £ u tr e 

1 

T 



- £ u e ii - ^(73^ + 72 tre) £ u 6> 



-1/T)KW 
+ fe" V"(-/x/T) 



(— 2Xi/T) £ u £( t iu) 
(q) 

-X2/T) £ u e^ 



+ [ire 6 (-1/T) trif 



(2.33) 



a/) . 

\(1/D) (trr-trr (q) ) / 

Thus, if we require that each term be separately positive definite, we obtain the following 
equations: 



-(2X 1 /T) £ u e^ u ) 
■(A 2 /T) £ u e^ 



2 {Q + v) 



tn + tr) 



( 



-ami 



(K + C) 



-VT)K {M 

h;d„(-n/T) 

( tT£ \ 
e 

{(-1/T) tvKj 



£ u (tie)\^ 
£ u 6 J 
V (V^)(trr-trr (q) ) / 
Here Q , H , and /C are antisymmetric matrices, 

/ £' /C \ 

-K! -/Ca 
-/C /Ca 

and r\ , cr , and £ are positive semidefinite symmetric matrices, 



(2.34) 
(2.35) 

(2.36) 



Q 
-G 



n 
-u 



7C 



(2.37) 
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»7i ??2 

J?2 773, 



cti a 2 

,0"2 0"3, 







C2 


cA 


c = 


C2 


Ca 


c 5 




u 


Cr, 





(2.38) 
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Note that only the symmetric matrices contribute when substituted to the entropy produc- 
tion rate f !2.33p . This means that the matrices rj , a , and £ are associated with irreversible 
processes, while the matrices Q , H , and K are with reversible ones. 

The relationship between the equations given above and the corresponding ones given in 
[T5] is summarized in Appendix |A] 



2.3. Fundamental equations 

Using Eqs. (I2.34p -f l2.38p at each point x = (x° = t, x) on timeslice S t , we can express (A) 
the currents and and (B) the evolution of strains, £ u e^ v ) , £ u ^^ , <£ u tre , and £ u 9 , 
only in terms of local thermodynamic quantities on E t . 

We thus conclude that the dynamics of relativistic viscoelastic materials is described by 
the following two sets of equations [131 IZ] : 

(A) current conservation laws: 

V M T^ = V M (e uV + t^) = , (2.39) 

= V M (n + z/ M ) = , (2.40) 

with the constitutive equations 

{K - Q tie- {Ka + Cs) 9 + ^ tr , (2.41) 



v*= -(n-a 2 )e^ + a 3 h^d u ( K -^y (2.42) 
(B) rheology equations: 

* , (2.43) 

, V^(-f), (2.44) 




-CiTtr£-(/C' + C 2 )T# + (/C + C4) tiK 
(£' -( 2 )Ttre-( 3 T6-(lCa-( 5 ) trK 

(73C4+72(/C'-C 2 ))T tj . (72C3-73(/C ; +C2))r q j 73(/C+C4)+72(^Ca-C 5 ) tr ^ N 

dct 7 det 7 dct 7 

(72Ci+7i(/C , -C2))T t _ (7iC3-72(/C'+C2))T « _ 72 (/C+C4)+7i (Ka-Cs) t ^ 
dct 7 dct 7 det 7 



(2.45) 
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The former set of equations describes the dynamics of D + 2 conserved quantities (p^ = 

eu^, n) , while the latter that of D(D + l)/2 dynamical variables = (£^,£^,0) . 

It is convenient to introduce the following parameters: 

r s =^, t, = (2.46) 

2det7 , . 

r± = — , (2.47) 

T (P C7 T ^ - 4 det 7 (det C s + /C' 2 )) 

a ± s -2«3^- (^ + C a )^) > (2 _ 48) 



Ca 7i - Ci 73 - 2/C' 72 ± ^ - 4 det 7 (det C s + /C' 2 ) 

where = (371 + Ci73 ~~ 2C272 > , and £ s is the principal submatrix of £ defined by 
C s = ( ^ § ) " Since C i s positive semidefinite, det £ s is non-negative. Note that r s , r a , and 
Rer-t are all non-negative. We further introduce the scalar variables 

e ± = i(tre-a ± 0). (2.49) 

Then the rheology equations f l2.43l) - fl2.45l) can be rewritten in a more compact form: 
(B') rheology equations: 

£ u £{^u) = £ (fiu) + —7—^- K (tiu) , (2.50) 



1 (ft + <r 2 ) T 



TVr A 



2 



« 1 e , (/Ca - Cs) (g± 7i + 72) + (/C + Q (a± 72 + Ts) , ~ ,„ , x 

r ± 2 det 7 

From these, we see that r s , r CT , and Re r± give the typical time scales for the relaxation of 

strains. 

The relation between the viscoelastic models and a few well-known rheological models 
(such as the Maxwell model and the Kelvin- Voigt model) is discussed in Appendix [B] 



3. Fluid and elastic limits 

In this section, we discuss the limits of elasticity and fluidity in the relativistic theory of 
viscoelasticity. We first identify the properties that characterize a given material as a fluid or 
as an elastic material. We then consider the long-time and short-time limits of our dynamical 
equations and show that fluidity is universally realized in the long time limit. We also make 
a comment on the subtlety existing in Maxwell's definition of viscoelasticity. 
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3.1. Fluidity and elasticity 

Fluidity is characterized by the property that the relaxation of the strains E^ u = (e^, e^, 6) 
proceeds instantaneously. Thus, their rheology equations are expressed as 

£ u £^ = , £ u e fl = , £ u 6 = , (fluids) (3.1) 

or equivalently, 

£ u £(nv) = , £ u s fl = , £ u e± = . (fluids) (3.2) 

This situation can also be realized in the long time limit, and we show in the next section 
that the constitutive equations for our viscoelastic model universally reduces to those for 
the Navier-Stokes fluids in the long time limit. 

On the other hand, elastic materials by definition do not undergo any plastic deforma- 
tions, and thus their intrinsic metric does not evolve for any processes. Thus, a given 
viscoelastic material is regarded as being elastic when its rheology equations are expressed 
as [13 [7] 

= - £ v h iLV = . (elastics) (3.3) 



3.2. Long time limit as a fluid limit 



Let the time scale of observation be T Q b s • If the observation is made much longer than the 
relaxation times (i.e., T Q b s S> r s , r a , Rer±) , then we can neglect the terms £ u £(^u) , £ u ^n , 
and £ u e± in Eqs. fl2.50p -f l2.52p because, for example, £ u e^ u ) ~ T~^ s e^ v ) <C r" 1 eu v ) . We 
thus obtain 



Q + ??2 
Al 



K 



Q + m 



K 



e± ~ t± 



2 v ± , a l \ T 

(Ka - ( 5 ) (a± 7i + 72) + (£ + C4) {a± I2 + 73) 



(3.4) 
(3.5) 



trK 



tie 



or 



2det7 

' C3(C4 + /C)-(C 2 + /C')(C5-M 



(3.6) 



J (detc s + /c' 2 )T \^_ (C2 _x:')(C4 + x:) + Ci(C5-x:a) / 

By substituting these equations to Eqs. (12.41 j) and f )2.42p . the constitutive equations take 
the following form: 

^(long) = r (q) _ K _ (tr K) ^ ? (3J) 
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v, 



(long) 



where we have defined viscosity and diffusion coefficients by 
det rj + g 2 



Cns 



rjiT 

det C + K 2 (a 2 Ci + 2a C 2 + Cs) - 2/C /C'(aC 4 + Cs) + K' 2 Ce 



(det C s + /C' 2 ) T 



det a + H 2 



(3.8) 

(3.9) 
(3.10) 
(3.11) 



Note that they are always non-negative, as can be seen from the inequality 

//Ca\ 

K? (a 2 Ci + 2a( 2 + Cs) - 2/C /C'K 4 + Cs) + £' 2 Ce = (/Ca K, -K 1 ) C K. 

\-K>) 



> 0. (3.12) 



In particular, when the material is locally isotropic, we can take r^K = P h^ u , with P the 
pressure, and thus the stress tensor certainly gives the constitutive equations for a relativistic 
Navier-Stokes fluid: 



_fJ,U 



, (long) - -2 Vns + (P - Cns tr K) hT . (3.13) 
We thus confirm that our viscoelastic model always exhibits fluidity in the long time limit. 



3.3. Short time limit as an elastic limit 



In contrast, at short time scales (T obs <C r s , Re r±) , we have 



1 1 

£ u £( f iu)'^ > £(uu) > £ u £± > e±j (3-14) 

so that Eqs. f)2.50p -( 12.52p can be approximated as 

£ u £{v,v) — ~ v ^ 2 ^< M ^> ; (3.15) 

^ (^-C 5 )(q ± 7i + 7 2 ) + (/C + C4)(q ± 7 2 + 73) (3 16) 

± 2det7 ' V 7 

> A(tre) * (^-C5)72 + (/C + C4)73 ^ (317) 

det 7 / 

By substituting Eqs. fl3.15p - fl3.17p into Eq. f l2.4ip . the stress tensor can be rewritten in the 
following form: 

r (short) = (q) _ 2 (G- m ) £/ X £ El \ 

' fiv r iw I 12 }t (jus) ^_|_^jj X,lfc W 
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(/C-CO tre- (/Ca + Cs)^ 

+ T7 Ce_det7 j 

[(/Ca - Cb) 72 + (/C + C 4 ) 73] T ul ; J MI/ 1 J 

These constitutive equations have the same form as those of a Kelvin- Voigt material (see 
Appendix |B]). However, one cannot yet identify the material at short time scales with a 
Kelvin- Voigt material, because they generically obey a different type of rheology equations. 

As discussed in the first subsection, elasticity is characterized by the condition that the 
intrinsic metric h^ u does not evolve, and the rheology equations for elastic materials are 
given by = , or equivalently by £ u £/j,v = K^ u [61 [151 E]- However, this is realized only 
when the conditions Q + 772 = Ai and (/Ca — (5) 72 + (JC + (4) 73 = det 7 are satisfied. That 
is, for generic values of parameters, even if the observation time is sufficiently shorter than 
the relaxation times, the intrinsic metric evolves when the induced metric h^ y does (i.e., 
Kfiv 7^ if K^y 7^ 0). Thus, Maxwell's original definition of viscoelasticity (considered only 
for the situations where the induced metric is static, K^ v = (1/2) £ u h liv = 0) needs to be 
modified for generic values of parameters, such that h^ u is allowed to evolve when does. 



4. Simplified Israel- Stewart fluids 

In this section, as an interesting example, we consider the case where JC' = 773 = 03 = £g = 
and = P h^ u . In this case, from the positivity of matrices r) , cr , and £ , the conditions 
772 = 02 = C4 — Cs = a l so must be imposed. Then the conserved currents take the following 
formal 

= e u"u v + Ph> w - 2Q e { " v) - JC (tr e - a 9) h" u , (4.1) 
n» = nu» -He^ (4.2) 

and the rheology equations become 
1 Q 

£ u £(nv) = + t K(jju/) 1 (4.3) 

£ u e, = --e, -^h»dJ-t) , (4.4) 



7"<7 A 2 " "V T 

Utre) = - (73Cl ~ 72C2)r tre + ^~^ T 9 + ^ ( ° 72 + 7s) tr K , (4.5) 
det 7 det 7 det 7 



6 From this form of the bulk stress and the relation 8 ~ (T — T)/T , we see that a/T can be identified 
with the thermal expansion coefficient. 
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£ u 9 = (72Cl ~ 7lC2)T tre - (7lC3 ' 7aC2)T 9 - /C( 7 1 + 72) tr K . (4.6) 
det 7 det 7 det 7 

By using the relations 

r (H = -20^), ^ = -Ue», n=-(trr-trr (q) ) = -X:(tre-a0) J (4.7) 

the rheology equations can be rewritten as 
1 2£ 2 

At (/if) = 7" ( At! /) 7— , (4.8) 

^ = - ^ + ^ V ^ (-£) > ( 4 - 9 ) 

£ n = _ [(« 72 + 73) Ci ~ (a 7i + 72) C2] T JC 2 (a 2 ^ + 2a7 2 + 73) 

det 7 det 7 

/CT[aCi(a7 2 + 7 3 )-C 2 (a 2 7i-73)-C3(a7i + 72)] . 1fV . 

+ d^ (4 - 10) 

£„9 = (7lC2 ~ 72Cl)T n - ^ Q7l + 72 ) tr X 



/C det 7 det 7 

[71 (a C2 + Cs) - 72 (a Ci + C2)] T 



9. (4.11) 



det 7 

This model gives hyperbolic differential equations for small perturbations around a hydro- 
static equilibrium, as is shown in Sec. [5j 

For brevity, we here consider the case when 9 is decoupled from other variables. This can 
be realized by setting a = 72 = C2 = in the above equations, and the rheology equations 
become 

t s £ u T(fiu) = - r^ v ) - n NS K^ u ) , (4.12) 

r a £ u v lx = - + cr NS hp d v (-^J , (4-13) 

r b £ U U = -U-C NS tTK, (4.14) 

£ u 9= -^9. (4.15) 

73 

Here we have introduced r b = 7i/(Ci T) , and the viscosity and diffusion coefficients are 
given in this case by r] NS = t s Q 2 /X 1 = Q 2 j{r]{T), ( NS = r b /C 2 /7i = K- 2 /((iT), and 
a NS = "H 2 /o"i. These equations look like the nonlinear causal dissipative hydrodynamics 
proposed in [13]. Although the nonlinear terms in [T4] (e.g., h^p v V p u v ) are important 
for numerical simulations of ultra-relativistic dynamics, these terms, in principle, cannot be 
treated properly in our first-order formalism. However, if we do not make the approximation 
£ u (txe) ~ tr(£ u £ Miy ) , then Eq. (I4.14p becomes — r b /C tr(£ u e fiu ) = r b (£ u Il+(l/D) trfTn — 
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JCK^ U ' e (^iu)) = — n — Cns trif and coincides with Eq. (14) in [2] where the spatial dimen- 
sion is set to be D = 1 . 

If we neglect the nonlinear terms, we then get relations of Maxwell-Cattaneo type: 

= _ 2 Vns - Ts h vS u p V p tt 7(5 , 

n = -Cns trK-r b n 7 V 7 n, (4.16) 
^ = a NS W d v (- 1 ) - r a h\ u"< VX , 

where ~k^ v = — . They are the constitutive equations for the simplified version of 
the Israel-Stewart modelJll 

Thus, in this case the rheology equations are equivalent to the constitutive equations 
for the simplified Israel-Stewart model P~TB"]) . and the [D + 1 + 1 + D(D + l)/2 + D] 
dynamical variables (excluding 9) can be determined from the D + 2 conservation laws 
(V^ = V M T^ = 0) and the D{D + l)/2 + D equations (OBI) . 



5. Hyperbolicity and dispersion relations 

In this section, we study linear perturbations around a hydrostatic equilibrium in Minkowski 
spacetime. We exclusively take a coordinate system (x M ) = (x°, x l ) in which the background 
metric is written as g^ u = r]^ u = diag(— 1, 1, ■ • • , 1) . A hydrostatic equilibrium is then 
specified by the velocity U( ) = u^d^ = d (i.e., vJL^ = 5$), the proper energy density e( ) , 
the number density rim) , and the vanishing strain tensor E$ = . The induced metric is 
then given by h$ = rj^ + uffluv = diag(0, 1, ■ ■ • , 1) . Note that from the fundamental 
relation for the hydrostatic equilibrium, §( ) = <7(o) (e(o) , W(o) , v/fyo)) = y/h{0) s (o)( e (o); n (o)) > 
other thermodynamic quantities such as the temperature Tm) , the chemical potential fim) 
and the pressure P(o) are determined as 

53(p) = J-tfe (0 ) - ^r^(o) + ^p-sJh^, (5.1) 
1(0) l(o) 1(0) V 



or 



5s (o) = 5e (o) - ^ 5ri( ) (5.2) 

7 (0) ^(0) 



7 The constitutive equations for a simplified Israel-Stewart fluid is obtained by setting the viscous-heat 
coupling coefficients to be zero in those for an Israel-Stewart fluid (i.e., ao = ai = in Eqs. (8a)-(8c) in 
ED- 



IT 



with the Euler-Gibbs-Duhem relation 

'(0) /^(O) , P(0) 



'(0) 



T(p) T (0 ) T (0 ) 



(5.3) 



5.1. Linear perturbations around a hydrostatic equilibrium 

We now consider linear perturbations around the hydrostatic equilibrium, 

9txv = Vfiv + > uil = + ^ ; V' = ^ + Von 8u u + 5u M 7/o„ , (5.4) 

e = e (0 ) + Se , n = n (0 ) + Sn , £ M „ = + E^ u , (5.5) 

and denote their conjugate thermodynamic variables by 

T = T {0) +5T, (i = m Q) + 6n, P = P (0) +5P. (5.6) 

We only consider the locally isotropic case: = P h^ u . Using the identity — 1 = = 
— 1 + 25uq = — 1 — 25u° , we can show that 5uq = 5u° = 0, and the acceleration vector 
a M _ u v Q vU v> — d Q $u^" has only spatial components: a = OqSu = and a 1 = do5u l . 
Moreover, from = e^ v u v = s^q , also has only spatial components, , in this linear 
approximation. Similarly, since = K^ v u u = K^q , the extrinsic curvature also has only 
spatial components, which are expressed as 

1 1 

K ij = 2 hitf (d^u u + d u u^) = - (di5uj + diSuj) , (5.7) 

or 

tr K = <W , K {ij) = l - (diSuj + djSui - ^ (d k 5u k ) . (5.8) 

As for the stress tensor (12.411) . by decomposing it as t^ v = rj^J + 8r^ v , the zeroth part is 
given by = P( ) , and from = r M „ u u = 5u l + <5r M o we can show that 5tqo = , 
5r i0 = —rfj 5ui = —P(o) Su,i , and the spatial components are written as 



(c 

2 



Sry = 5P -2(G- m ) e {ij) - ^- faduj + djSm - ± (d k 5u k ) h m 

i(o) L ^ 



- [(£ - C 4 ) tr e - {Ka + Cs) 0] hf - A (fl^*) fcg> . (5.9) 

J (o) 

The diffusion current is written as 

^ = -{U- a 2 ) e" + a 3 <M(-|) . (5.10) 
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We now substitute the above expressions to the set of fundamental equations, consisting 
of (A) the conservation laws (12731?]) - ([2732]) and (B) the rheology equations (12741?]) - ([2735]) (or 

(A) As for the conservation laws of energy-momentum tensor, the component along is 
given by = u v d^T^ lu = d^(T^ lu u u ) — T^ LV d^u u = —d^{eu^) — t^ v d^u v . From this we obtain 

<9 M (eu M ) = d 5e + e( ) di8u l 

= -i*" d^u u = -7$ d^u v = -P (0) dM , (5.11) 



or 



d Se = —iO(o) di5u l . (5-12) 

Here = e(o) +P(o) is the enthalpy density at the hydrostatic equilibrium. As for the com- 
ponents orthogonal to , from the equations = h\ v d^T^ u = h\ v d^eu^v?) + h\ v = 
e d^ux + h Xu d fM r^ v = e a\ + h x v d^r^ , we obtain 

e di = e( ) d 5ui 

= -K 9%, = -VSt^ = -d 5r 0l - d k 5r lk = P (0) d°5 Ul - 8 k 5n k , (5.13) 



or 



W(q) d 6ui = - d k 5T ik 

■(£>- 2)7/3 , Ce \ V3 



- d^P + 2(Q- V2 ) d k e m + ( [ nrp )V3 + i?-) <Wu fc + A<W 

v ^ J (o) J (o) y J (0) 

+ (/C - C 4 ) 9 t (tr e) - {Ka + Cs) $0 , (5.14) 



where A is the spatial Laplacian, A = <5 U di dj . The conservation law of particle number 
current becomes 

= a M (nu" + z/ 1 ) = d 6n + n (0) <W + <r 3 A <*(-^) - (W - <r 2 ) . (5.15) 
(B) The rheology equations are linearized as 

doS(ij) = ^ {p i5u i + d 3 Su * ~~ ~]j ( d k5u k ) h[f^j - ^e {ij) , (5.16) 

_i £i _(*±p)Zk Sj *(_a ( 5.i7) 

a v (73 Ci + 72^C' -C2)) T (0 ) (72 Ca - 73 (£' + C2)) T (0 ) 

9 (tre) = — treH — * 9 

det 7 det 7 

19 



7 2 (JCa - C5) + 73 (jC + C4) 
det 7 



(5.18) 



« fl (7i (K' - C2) + 72 Ci) T (0) (71C3-72 (g + C2)) T (0) 
= — — — tr e r- — — 



det 7 

72 (JC + (4) + 7i ~ (5) 
det 7 



det 7 



(5.19) 



where we have used the approximation £ u e^ ~ d^e^ , £ u Ei ~ <9 £; , £ u {tie) ~ 9 (tre) , and 
£ u ~ . 

Since we are considering locally isotropic materials, the fundamental thermodynamic 
relation ( I2.24p can be rewritten with the use of the Euler relation (1C.6|) as 



5s 



Se 



5n-- 2Ai e<""> 5e^ u) - - ( Tl tr e + 72 0) 5(tr e) 



T 
1 



A 2 £"fe M -- (73^ + 72 tre) <W. 



(5.20) 



If we denote the thermodynamic variables collectively by (a r ) = (e, n, E( uu ), e n , tr e, 6*) , the 
matrix A = — {d?s/da r da s ) L Q , is positive definite from the convexity of entropy. Here |(o) 
means that the matrix is evaluated at the hydrostatic state. In the following discussions, we 
assume for brevity that the matrix takes the following form: 



A 





A 2 











o\ 


A 2 


A 3 




























































A 6 


A 7 


v° 











A 7 


a 8 / 



(5.21) 



where the principal submatrix 



A. 



(A, 


A 2 \ 


\A 2 





/ 


d 2 S 




d 2 s 


) 




( 




de 2 


(0) 


dedn 


(0) 








d 2 s 




d 2 s 








V 


dedn 


(0) 


On 2 


(0) ) 




V 



de 
d(l/T) 



dn 



(0) 
(0) 



d(M/T) 
de 

dn 



\ 

(0) 
(0)/ 



(5.22) 



is positive definite. Then the Gibbs-Duhem equation ( 1C.7P can be written ag_ 

diSP = s( ) di5T + n (0 ) 

= S(0) d^iii/ry 1 ] + n (0) 9, 5 [(l/T)- 1 o*/r)] 



8 Note that the right-hand side of Eq. (|C.7|) can be set to zero for the linear perturbations around a 
hydrostatic equilibrium. 
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(w(o) Ai + n (0) A 2 ) 2(o) <9;5e + (u> (0 ) A 2 + n (0 ) A 3 ) T (0 ) ft<5n , (5.23) 



di6(-£) =-(A 2 d i 6e + A 3 a i 6n), (5.24) 
and we finally obtain the following set of linearized equations of motion: 

d 5e = — W(p) di5u l , (5.25) 

w {0) d 5 Ul = 2{Q- Tfc) d k e m + f (D ~ 2)7/3 + W«* + ^- Atfu* 

v i ^ i (o) J (o) y J (0) 

+ (X:-C4)9i(tre)-(/Ca + C 8 )^ 

- (w( ) Ai + n( ) A 2 ) T( ) <9;<5e - (u>( ) A 2 + ri( ) A 3 ) T (0 ) <9j<5n , (5.26) 
d Sn = - n (0) S^u* + a 3 (A 2 A<5e + A 3 A5n) + (ft - <r 2 ) fl^e* , (5.27) 

d £(ij) = ^ { diSu i + d i Sui ~ J) W uk ) h if) ~ \ £ m » ( 5 - 28 ) 
(fe = - - Si + ( ^ + ^ )T(0) (A 2 dje + A 3 c^n) , (5.29) 

T a A 2 

fl s (73 Ci + 72 (£' - C2)) ^o) (72 Cs - 7s (/C' + C2)) 2 1 (0) 

<9 (tr £j = - — tr e + - — - '- 9 

det 7 det 7 

+ 72(/Ca-C 5 )+73(X: + C4) ^ s (5>30) 
det 7 

a a ( 7l ~~ + 72 ^) T (°) * ( 7l ^ 3 ~ 72 ( /C ' + &)) T (°) a 

d 6 = ± — '- tre - — - 6 

det 7 det 7 

- 72(/C + C4) + 7l(/Ca ~ C5) ^. (5.31) 
det 7 

We now consider wave propagations in the x D direction, demanding that perturbations 
depend only on x° and x D : 

5 Ul = 6 Ui (x°, x D ) , e i:j = Eijix , x D ) , 5e = 5e(x°, x D ) , 5n = 6n(x , x D ) . (5.32) 

Then the above equations can be rewritten as follows: 

d £(ii) = - ~j^t~~ d D 5u D - — e^u) , (5.33) 

JJAi T s 

d e {IJ) = £(/j) (I^J), (5.34) 

doe (id) = ~ ^J 12 d D 5uj - — e {ID) , (5.35) 

iO(o) doSuj = 2(S- m) d D £u D ) + <P D 5ui , (5.36) 

J(0) 

= Ei , (5.37) 
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W(o) d 5u D = 2(Q- f] 2 ) d D e {DD) + ( 2 ^ D + d 2 D 5u D 

\ L) i rm J cm / 



<9o<5e = —W(q)8d5ud, (5.38) 

-2(£>- 1)773 Ce 

n ' rri 

(0) 1(0)- 

+ (/C - U) 9u(tr e) - (/Ca + Cs) <9d# 

- (w(o) Ai + n( ) A 2 ) T(o)d D 5e - (i0( ) A 2 + ri( ) A 3 ) T^d D 5n , (5.39) 
<9 5n = - ri( ) d D 5u D + cr 3 (A 2 d^tfe + A 3 d 2 D 5n) + (H - a 2 ) d D e D , (5.40) 

(£)-l)(g + 7 fc) l 

<9 £(dd) = prr o D du D e {DD) , (5.41) 

L> Ai t s 

<9 £d = - - e D + (H + f 2)T(0) ( A 2 «9^e + A 3 d D 5n) , (5.42) 

Tcr A 2 

s (73 Ci + 72 (£' - Ca)) 7(o) (72 Ca - 7s (£' + C2)) T (0 ) 

do(tre) = — tre + ^ — 9 

det 7 det 7 

72 (£a - Cs) + 73 (£ + C4) a , 

H r- cMwd , (5.43) 

det 7 

„ a in (r - c 2 ) + 72 Ci) r (0 ) , (71 C 3 - 72 (/r + c 2 )) t (0) 

a 6> = — — tie — — — 

det 7 det 7 

72 (jC + C4) + 7i (/Co - Cs) a r ,k„„n 

— o D bu D , (5.44) 

det 7 

where /, J = 1, ■ ■ ■ , D — 1 . This set of equations can be further decomposed according to 
the transformation properties under the little group SO(D — 1) : 

1. tensor modes: (e/jn, ; 

2. shear modes: (£{id), Sui, £j) ; 

3. sound modes: (tre, epu), 5m.d, 5e, 5n, £d, 0) . 

In the remainder of this section, we study hyperbolicity and dispersion relations for each 
type of perturbation modes. 



5.2. Tensor modes 

For tensor modes, the set of equations can be written as 

doe (ii) = ~rT\~~~ d D8u D - — £<//) , (5.45) 

JJAi T s 

d £(u) = £(U)- (5.46) 

7~s 
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From the identity ^ 7 e^i) +£{dd) = , the number of independent variables of E(n) is D — 2 . 
If we define the variables E u by 

{£(//) — e ' Ud-\)(d-\)) (for I = J) 
(5.47) 
e</j> (for / ^ J) 

then the number of independent Eu is D(D — l)/2 — 1 = (-D — 2)(Z) + l)/2 because 
E(d-i)(d-i) — 0, and the equations for Ejj become 

= d E u + -Ejj . (5.48) 

Thus, if we consider plane waves propagating in the x D direction, 

5u t = &Ui{uj, k) e lkxD - iulx0 , Ei j = eij{uj, k) e lkxD - iux ° , (5.49) 

we obtain the dispersion relation u = — i/r s which represents nonpropagating, purely dissi- 
pating modes. Since r s is positive, the imaginary part of u is always negative, and thus we 
find that the tensor modes are always stable. Such relaxation modes correspond to stress re- 
laxations observed at rheological time scales (T obs ~ r s ), and will disappear at hydrodynamic 
time scales (T obs ^> r s ). 

5.3. Shear modes 

For shear modes, we have the equations 

= (do + T' 1 ) e {ID) - d D 5uj , (5.50) 

= <9 5u/ - — — — d D en D \ d 2 D bui , (5.51) 

W( ) W(0) (o) 

0= (do + r^ej. (5.52) 

Note that Ej is decoupled from the other variables, and Eq. ( 15.521) represents its pure relax- 
ation with relaxation time r a (> 0) . 

If we set 7/3 = , by redefining the variables by 



si± = ±\I—^-£(id) + t;8u i , (5.53) 
iU(o) 2 



the set of linearized equations for si± can be written as 



d =F Qhcar d D ) s I± ± Sl+ 2r Sj = , (5.54) 
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for I = 1, 2, . . . , D — 2 . These are hyperbolic equations and the characteristic velocity is 
given by 



Cshear — \ l ■ (5.55) 

For generic cases, from Eqs. ( I5.50p and ( I5.5ip . we obtain telegrapher's equations with 
Kelvin- Voigt damping 



(d 2 + -d - did, - c s 2 hcar 91) ( £{ID) ) = . (5.56) 

V r s W( ) T( ) / \ Sm I 



Although they are generically nonhyperbolic and have infinite wave-front velocity as in 
the standard relativistic fluid mechanics, they can be made into hyperbolic telegrapher's 
equations by setting rj 3 = E 

If we consider the short time limit (r s — > oo), the differential equations become 



W(o)l(o) > \5ui 

The wave equations in this form also appear for viscous solids such as Kelvin- Voigt materials, 
and reduce to the wave equations when r] 3 = . 

Finally, for plane waves 

6ui = ^{lo, k) e lkxD -'^ x ° , Ei j = e^uj, k) e ikxD -^ x ° , (5.58) 
from (I5.56p . we obtain the dispersion relation 



T 2 + - + — ^- k 2 )T + c hear k 2 = , (5.59) 

where T = —iu . Since all the coefficients are positive, the real part of T (or the imaginary 
part of uj) always takes negative values, and thus we see that there are no unstable growing 
modes in the shear modes. Equation (I5.59P has two solutions, which are expanded around 
k = as 

u = | "i + i (1 - r.) dear ^ 2 + i (1 - r.) 4 ear r s 3 k* + 0(k«) 

_ i C shcar T s k 2 — 1 (1 — T s ) C 4 hear T s 3 k 4 + 0(k 6 ) , 



9 In this case, from the non-negativeness of the matrix 77 , ry 2 (and thus det r/) must vanish. However, this 
still gives a positive shear viscosity if Q ^ , as can be seen from Eq. 
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with r s = 773/ (rj NS T^) . The former represents the relaxation modes which are not observed 
at hydrodynamic time scales (T obs 3> r s ). The latter represents the hydrodynamic modes 
where u — > in the limit k 2 — > , and from the coefficients of k 2 , the diffusion coefficient is 
found to be c 2 hcar r s = rj NS /w( ) . Moreover, by the comparison with the dispersion relation 
of Maxwell-Cattaneo type, the effective relaxation time associated with the hydrodynamic 
modes is read off from the coefficients of k A as (1 — r s ) r s . Indeed, if we set r s = 1 , the 
effective relaxation time becomes zero and the dispersion relation becomes purely diffusive; 
u = -i (rj NS /w { o)) k 2 . 

If we are interested only in the hydrodynamic modes, the dispersion relation coincides 
with that of the Israel-Stewart model up to 0(k 4 ) by identifying (1— r s ) r s with the relaxation 
time r n in the Israel-Stewart model. However, if the relaxation modes are also taken into 
account, our viscoelastic model has a richer structure than the Israel- Stewart model, which 
is the special case (r s = 0) of the viscoelastic model. 

5.4. Sound modes 

Finally, for sound modes, we have the following set of differential equations: 

= d 5e + W(p) d D 5u D , (5.61) 

= d 5u D d D e {DD) - — — - + — d D 5u D 

w {0) V D w (0) T (0) iu ( o) T (0) / 

o D (tr e) H d D 9 

W( ) W( ) 

+ K 0) A 1 + n (0) A 2 )T (0) ^ + K)A 2 + n (0) A 3 )T (0) ^ 

= d Q 5n + n {0 ) d D 5u D - a 3 (A 2 d 2 D 5e + A 3 d 2 D 5n) - (H - a 2 ) d D e D , (5.63) 

(D-l)(g + T)2 ) 1 
= OqS (dd) prr o D 5u D + — e {DD) , (5.64) 

= d e D + —s D - & + P- T(0) (A 2 d D 5e + A 3 d D 5n) , (5.65) 

n a u x , (73 Ci + 72 (£' - C2)) T {0) 

= a (tr£:) + — — trs 

det 7 

72 (JCa - Cs) + 73 (JC + 0) a . (72C3-73(/C / + C2))T (0) /a 
del^ ^ det^ (5 - 66) 

n Qfl 1 72(/C + C4)+7i(/Ca-C 5 ) a , 
det 7 
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Mi = 



M A = 



<J C 2 (072 + 73) 
W(o) det 7 

H (w {0) A 2 + n (0 ) A 3 ) 



Mo 



/2 (£> - 1) ?7 NS 



D W(q) r s 



M, 



/ / C 2 o (071 + 72' 
it)(o) det 7 



(0) 



/ 



M~ l B x M 



(73 Cl+72 (^ / -C 2 ))T(Q) 



dct -7 




" 






((K / -C 2 )7l+Cl 72)^(0) 



Mr = — 



00000 

i 

o s 
00000 
00000 
0000 — 



•H / det A s T( ) u>( ) 



(5.76) 



(CC / +C 2 )73-C3 72) r (Q) / Q71+72 " \ 



dct ~y 



a(a72 +73) 










l(a72 + 73) 

a 71+72 



(71 C 3 -72 (K'+C 2 ))T (0) 



dct ~y 



J 



(5.77) 



The real matrix M 1 B Q M is symmetric and can be diagonalized. The eigenvalues are 
calculated to be {0, 0, 0, ±v±} , where 



„2 I Cns I 2(D-1)t? N s I A 3 CT NS 

^ D 1U (0) T S f T(j 



"(0) T> 



'-\i(c 7 \ CnS I 2(D ~ 1),7 NS I A 3 g NS V 
! VV S W (0) T b Dw (0) Ts T " J 



4A 3 CT NS [ CNS 1 2(P-1)»7NS 



™(0) T b 



dct A s m(Q) T (Q) 



(5.78) 



give the characteristic velocities. Since all the eigenvalues are real, we see that the system 
of differential equations f !5.68p is hyperbolic. 

If we particularly set H = (and thus <r NS = 0), the characteristic velocity reduces to 



v± 




2(£>-l)7fc 



W( ) r b 



Dw 



(0)T B 



(5.79) 



and agrees with the large wave-number limit of the group velocity (which in our case coincides 
with the front velocity and the characteristic velocity) in the Muller-Israel-Stewart theory 
(see, e.g., Eq. (49) in [IS]). If we take the long time limit, Tb , r s — > , the characteristic 
velocity becomes infinitely large, and thus causality gets violated. 

For generic cases (i.e., when we do not impose the conditions 773 = (q = 03 — 0), from 
Eqs. f l5.6ip -f l5.67p . the dispersion relation for the plane wave 



5u t = 5u t (u,k) e ikx -™ x ° 
Se = Se(co,k) e lkxD -^ x ° . 



£ ii (cj,A;)e 1 ' 



ik x D —iujx 







Sn = Sn(u, k) e 



ikx D — iujx 



(5.80) 
(5.81) 



is obtained as 



r 7 + (c 60 + c 62 k 2 ) r 6 + (c 50 + c 52 k 2 + c 54 k A ) r 5 + (c 40 + c 42 k 2 + c 44 k A ) r 4 

+ (c 30 + c 32 k 2 + c 34 k 4 ) T 3 + (c 22 k 2 + c 24 k 4 ) T 2 + (c 12 /c 2 + c 14 A; 4 ) r + c 04 k 4 = , (5.82) 
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where r = — iu and 



c 60 =r~ 1 + T" 1 + + r_ 1 



C50 



C40 



C30 



C24 



C12 



0)4 



Ts r g + (r s +r g ) (t + +r_ ) +T + 7 
t s 7- ct t + t_ 



52 s ~ r ^(0)^(0) W(o) D TO(0) T (0) 

1 _l ^-1 1 ^-1^ 



+A 3 a 3 (r s - 1 + + r + x +r 
c 54 = A 3( t 3 (— Jl- + g=^p 



(O)-'(O) ^™(0) J (0) • 

+r g -+T + +r_ 



1 1 j Co 1 CNs(-r s -r CT + (r s +-r CT )r b 1 r + r_) 



C42=<(' 

4 



S \ S CT + — y tO( )T( )T s T CT ' U>(0)T B r CT T + 7 

2(D— l)r? N s(f B T CT +(f s +T CT )(r++T_)+T+T_) . A3 cr NS (t s f g + (r B +f CT ) (r + +t_ ) +r + t_ ) 



D«J(0)-r s T (T T + r_ T s T a T + T- 

r - A „ / ^(r^^r- 1 ) , CnS^ , 2(D-l)r ?3 (f B - 1 +f- 1 +r- 1 +rr 1 ) dctAT (0) ^ ( 
C 44 -A3^3l, T , n . Wn . + + Dr m «,„, + A3 



s (t b t ct + (r s +T CT ){t + +t_)+t + t-) Cns(t s +r CT +r b V+T- 



TsTo-T+T. 

£ _ c s( T s r g- + ^ | s~r'g-;^' + ~r , -;^ l + _|_ 
3 TsTctT+T- uij-ojTsTCTT-f-r- 

I 2(D-l)rj NS (f s +r J +T + +T-) A 3 q- ns (t b +f CT +t + +t_ ) 

= A CT ( Ce i_ CnS (t s t g + (r B +Tg )t~ 1 T+T- ) 

, 2(g-l)?? N s(f s f g + (f B +-r g )(T + +r-)+T + T-) 

-DlO( )T s T CT T-fT_ 

det AfT 3 T( 0) w (0 ) (f ~ 1 +r ~ 1 +r~ 1 +rj 1 
A 3 ctns 

<*> = ^l( c s(- S + r. + r + + r_) + ^ + + A 3 a NS ) , 



T S T CT T + T_ V » v " ' ' U!( ) ^™(0) 

A3CTNS f Cns (tb +f g +r~ 1 t + r_ ) 2(D-l)r? NS (f B +f tT +T + +T_) 

dct Aw( )T( 0) (r s f (T +(T s +f o -)(T + +T_)+r + T_ 



+ 



r s r CT r + r_ 



c _ A 3 g NS / Cns I 2(D-l)>y NS dct Aio (0) T (0) (r s +f g +T + +T_ ) x 



dct Aq- ns T (0) uj (0) 

T s To-T + T- 



Here we have defined non-negative constants 



??3 » . . t7 3 

^NSJ(O) CTnS 



and redefined Tb as 

Cns det 7 
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7+ = a 2 7i + 2a 72 + 73 > 0, (5.101) 
Pea = (Cs Ce - C5) 7i + 2 (C 4 C 5 - C2 Ce) I2 + (Ci Ce - C|) 7s > , (5.102) 

which becomes 7i/(Ci T) when the parameters are taken as in Sec. HI Note that complex 
parameters t± appear always through the combinations r + + t_ = 2Rer + (> 0) , t + t_ = 
\ r +\ 2 (> 0) or r+ + rZ = 2Rer + /|r + | 2 (> 0) . One can check that all the coefficients are 
positive, and thus at least the necessary condition for the stability is satisfied. For a full 
analysis to be performed, one should further check the Routh-Hurwitz stability criterion, 
which we have not carried out yet. 

The dispersion relation around k = gives seven solutions, and four of the seven take 
the following form: 

u = -— + 0{k 2 ), u = -- + 0{k 2 ), u=-— + 0{k 2 ). (5.103) 

They correspond to the relaxation modes, and as the observation time becomes much longer 
than the relaxation times Rer-t, r s , and r a , these modes fade away in time and will not be 
observed eventually. 

The remaining three modes are hydrodynamic modes and have the following expansion 
in k : 

r 2 \ 



u =c s \k\ -i Cl k 2 + (c 2 --J-) \k\ 6 + 0(k A ), (5.104) 

Z C s 



ci 



with 



Cl = i ( + 2 (-° ~ jj ^ + (A3r2( 0) + A 2 u>(o)) 2 a N5 T (0 ) - ^ 
1 2 \w(p) -Dw(o) c 2 w (0 ) 
c = / Cns (t+ + r_ - t^t+tJ) (D - 1)% S (1 ~ r s )r a 
\ 2to ( o) Dw{o) 

T (0) a NS (A 3 W(o) + A 2 W( )) 2 / 2 ( NS | 2(D-l)r) NS det A <r NS T (0) iU( ) 



2cl W(q) T a V S W(p) T a D W(p) T a C 2 T a 

(5.107) 

In particular, if we neglect particle diffusions ("H = <r NS = 0), we have 

ei = ^+ (D ~ 1)%S , (5.108) 

c 2 = f C - (r+ + r --^ V+T - ) + (D - 1) ^ s(1 - rs)rs ) c s . (5.109) 
V 2w(p) Dw {0 ) J 
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Up to 0(k A ) , this dispersion relation coincides with that of the Israel-Stewart model if we 
identify r + + t_ — r^" r+r_ and (1 — r s ) r s as the relaxation times Tn and r n of the Israel- 
Stewart model, respectively (see e.g., Eq. (47) in [15]) J 10 

6. Conclusion and discussions 

In this paper, we have studied the relativistic viscoelastic model [13] proposed recently on 
the basis of Onsager's linear regression theory on nonequilibrium thermodynamics. We first 
rederived the model using a local argument based on the current conservation laws and the 
positivity of entropy production rate. We then studied in detail the properties of the model 
and showed that our model universally reduces to the standard relativistic Navier-Stokes 
fluid mechanics if the observation time is sufficiently longer than the relaxation times. 

We also studied linear perturbations around a hydrostatic equilibrium in Minkowski 
spacetime. We showed that the wave equations for the propagation of disturbance become 
symmetric hyperbolic for some range of parameters, so that the model is free of acausality 
problems. This fact suggests that the relativistic viscoelastic model can be regarded as a 
causal completion of relativistic Navier-Stokes fluid mechanics, defining the latter as its long 
time limit. 

Although the wave equations are not hyperbolic for generic values of parameters, the 
problem of ill posedness in numerical simulations will be significantly remedied from the 
situations encountered in Navier-Stokes fluid mechanics. To see this, let us consider a shear 
mode as an example. As we saw in Sec. I5.3[ the dispersion relation in the long wavelength 
limit is given by Eq. f)5.60p . 

u = -i (r/ NS Mo)) k 2 - i (W^o)) 2 (1 - r„) r s A; 4 + 0(k 6 ) (6.1) 

and has the same structure as that of the Israel- Stewart model up to 0(k 6 ) so long as 
(0 <)r s < 1. This implies that, even for a parameter region where the wave equations 
are not hyperbolic, the behaviors at short wavelength scales are still remedied to an extent 
similar to that of the Israel-Stewart model, and thus the problems associated with the 
causality violation are expected to occur less likely in numerical simulations. It should be 
interesting to check this statement with a direct numerical simulation. 

10 In order for the correspondence to hold, we need to further choose the parameters such that r + + r_ — 
tC t+t- and (1 — r s )r s are both positive. 
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As discussed in Sec. El the dispersion relations for linear perturbations with generic 
parameters exhibit two kinds of branches. One is the "hydrodynamic branch," where u — > 
as k — > , and corresponds to the poles in retarded Green's function in the Kubo formula for 
dissipative fluid mechanics. If we neglect the effect of particle diffusion ("H = o~ NS = 0), these 
poles in the relativistic theory of viscoelasticity coincide with the poles of the Israel- Stewart 
model up to 0(k 6 ) for shear modes and 0(k 4 ) for sound modes [see Eqs. (I5.60p . (I5.106P and 
f l5.107p ] by identifying r + + r_ — t + t_ and (1 — r s ) r s with the relaxation times Tn and r w , 
respectively, in the Israel- Stewart model. In the so-called fluid/gravity correspondence [171 
T8] , such poles are actually found in retarded Green's functions calculated at the boundary of 
an asymptotically AdS geometry, and the relaxation time is obtained to have the value T n = 
(2 — In 2)/ (27rT) for strongly coupled M = 4 Super Yang-Mills theory. This suggests that we 
should set (1 — r s ) r s = (2 — In 2) / {2ttT) if we want to establish a mapping between the fluids 
described by strongly-coupled Yang-Mills theory and those described by our viscoelastic 
model. 

The other branch ( "rheological branch") gives a behavior that u converges to a non- 
vanishing, pure imaginary value, uj — > — i/r s + 0(k 2 ) , as k — > , and thus corresponds to 
the relaxation of strains. These relaxation poles are usually discarded in the discussion of 
viscous fluids, because the observation time for fluids is much longer than the relaxation 
times and the relaxation modes disappear at such time scales. However, if such poles can 
also be found in retarded Green's function at the boundary theory, then the fluid/gravity 
correspondence may be understood within a more general framework of the "viscoelastic- 
ity/gravity correspondence. "lj It would be interesting to pursue the study in this direction. 
It should also be interesting to investigate the viscoelasticity/gravity correspondence along 
the line of the recent study relating the solutions of the Navier-Stokes equations to those of 
the Einstein equations [191 ED] • 

As other future directions to be pursued, it should be important to extend the model such 
that one can treat more complicated systems like multicomponent viscoelastic materials. 
Such extension is actually straightforward and is under investigation. Another interesting 
direction is to extract the transport coefficients from kinetic theory or to extend the theory 
such as to include higher- derivative corrections. 



n To establish this, one first would need to investigate whether the parameters r s and r s can be obtained 
consistently for sound and shear modes. 
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A. Entropic formulation of relativistic viscoelastic fluid 
mechanics 

In this Appendix we give a brief review on how the fundamental equations [Eqs. ( I2.39p - 
(12.451) ] are obtained from the relativistic theory of viscoelasticity [13] constructed on the 
basis of Onsager's linear regression theory [9j [TOj [TTJ [12]. We use the same geometrical 
setup and the same definition of viscoelastic materials as those given in Sec. 12.11 See [TB] 
for a more detailed description. 

We assume that the local thermodynamic properties of the material particle at x (already 
in its local equilibrium) are specified by the set of local quantities (b A (x), c 1 (x), d p (x)^j . 
Here c 1 (x) denote the densities of the existing additive conserved quantities C 1 . b A (x) de- 
note the "intrinsic" intensive variables possessed by each material particle (such as strains), 
and d p (x) denote the remaining "external" intensive variables which further need to be in- 
troduced to characterize each subsystem thermodynamically (such as the background elec- 
tromagnetic or gravitational fields). We distinguish density quantities from other intensive 
quantities, and by multiplying them with the spatial volume element \/h , we construct new 
quantities which are spatial densities on each timeslice. For example, the entropy density s 
and the densities c 1 of conserved charges are density quantities, and for them we construct 
the following spatial densities: s = \/hs , c 1 = y/hc 1 . The local equilibrium hypothesis im- 
plies that the local entropy s(x) is already maximized at each spacetime point x and is given 
as a function of the above local variables; s(x) = s(b A (x), c 7 (x), d p (x)) . If we denote by 
(e s , e t ) the spacetime scale where the local equilibrium is realized, then at spacetime scales 
larger than (e s , e t ) , we need to take into account the effect that the material particles com- 
municate with each other by exchanging conserved quantities (such as energy-momentum 
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and particle number). The second law of thermodynamics tells us that, if boundary effects 
can be neglected, this should proceed such that the total entropy of the larger region gets 
increased. 

In order to describe such dynamics mathematically, we first introduce the spacetime scale 
(L s , L t ) which is much larger than the spacetime scale (e s , e t ) and assign to each spacetime 
point x = (x° = t,x) on timeslice £$ a spatial region E X [L S ] of linear size L s (see Fig. [3]). 
We then consider the total entropy of the region E X [L S ] : 




Figure 3: Time evolution of material particles in the large region S X [L S ] [T3] . 



S(t; £,[L S ]) = f d D ys(b A (t,y), c 7 (f,y), rf p (t,y)) . (A.l) 

The irreversible evolutions of intrinsic variables a r (x) = (b A (x), c 7 (x)) at x will proceed 
toward an equilibrium of the region E^Ls] . Due to the condition L s 3> e s , we can assume 
that the influence from the surroundings of the region S X [L S ] is not relevant to the dynamics 
of a r (x) because x is well inside the region. An equilibrium state of the region Ea.[L s ] will be 
realized when the observation is made for a long period of time, L t , and can be characterized 
by the condition 

Note that the functional derivative is taken only with respect to a spatial, D-dimensional 
unit in the functional. We denote the values of a r (x) at the equilibrium by a r Q (x\ L s ) = 
(b A (x; L s ), Cq(x; L s )) . One should note here that, since c 7 (t, y) are conserved quantities, 
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the variations (1A.2j) with respect to c/-type variables should be taken with total charges 
kept fixed at prescribed values: 

f d D yc I (t,y)=C / (E a; [L s ]). (A.3) 

J S x [L s ] 

A simple analysis using the Lagrange multipliers shows that the condition of global equilib- 
rium is expressed locally as 

ds 



A (x) = and h;(x)VMx)=0, (A.4) 
where j3j is the thermodynamic variable conjugate to c 1 that is defined by 

Pi(x) = -^(x) . (A.5) 

The total entropy of the region £ X [L S ] at an equilibrium is given by 

S (t; E X [L S }) = / d D y ~s{b£(t,y), <*(t,y), d p (t,y)) , (A.6) 

</£°[L s ] 

where S° [L s ] is a hypersurface orthogonal to the velocity field at the equilibrium, Uq ee Pq /e . 
When the material can be regarded as being at an equilibrium at spatial infinity, we can 
fix the labeling s of the new timeslices at the equilibrium with the labeling t of the 

original timeslices by setting s = t if E° conforms with H t at spatial infinity. If we 
denote coordinates corresponding to the new foliation by = (x'° = t, x' 1 ), then 
the velocity field uq = Uq d' will be expressed in the following form: 



This expression defines the new lapse N and the new shifts iVg at the equilibrium that 
are realized at spacetime scale (L s ,L t ). For configurations other than the equilibrium, the 
total entropy S(t; E X [L S ]) is smaller than that of the equilibrium So(t; E X [L S ]) , so that if we 
denote their difference by 

AS(t; E X [L S }) ee S(t; E x [L a ]) - S (t; S X [L B ]) , (A.8) 

AS* is always nonpositive. 

In the previous paper [13], it is proposed that the difference AS can be effectively written 



in the following form at the lowest order in the derivative expansion for linear nonequilibrium 
thermodynamics: 

AS(t; E X [L B ]) 
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Here the scalar function Nq is the lapse at the equilibrium defined in Eq. I I A. 71) . the coefficient 
{^^i /vtvjj is a symmetric, positive semidefinite matrix, and all the elements are spatial 
tensors, = = i^ I,uJ u u . The integral region can be replaced by £ X [X S ] because the 

difference is of higher orders in the derivative expansion. See the appendix in [T3| for a 
derivation of ( 1A.9j) for simple cases. The functional form of the total entropy, S(t; E X [L S ]) = 
£>o(t; E X [L S ]) + AS(t; H X [L S ]) , is called the entropy functional in [T3] . 

We now consider Onsager's linear regression theory [91 [101 [TTJ [12] assuming that the 
total entropy is given with this entropy functional. In Onsager's treatment the irreversible 
evolutions of thermodynamic variables a r (x) are given by 

[a r (x)] ilt = L rs f s (x). (A.10) 

Here f s (x) is the thermodynamic force defined by 

oa s {x) 

and in the relativistic nonlinear thermodynamics, the dot should be defined as a r = N £ u a r 
|13j . where £ u is the Lie derivative with respect to the velocity u = u fJi (x)d fJ- . V s are 
the so-called phenomenological coefficients and can be shown to satisfy Onsager's reciprocal 
relation EDI 03] 

L rs = (_ 1 )K|+|a s | L sr ^ (A. 12) 

where the index \a r \ expresses how the variables transform under time reversal, a r (x) — > 
(_i)k r l d r (x) o The Curie principle says that V s can be block diagonalized with respect 
to the transformation properties of the indices (r, s) under spatial rotations and the parity 
transformation [21], that is, under the subgroup 0(D) of the local Lorentz group 0(D, 1) 
in local inertial frames. For example, when a r constitute a contravariant vector, (a r ) = 
(a M ) , the equations of linear regression should be set for each of the normal and tangential 
components to the timeslice through x : 



5 AS 
_5a u (x)_ 



(A.13) 



12 When the background fields d p change as d p — > under time reversal, the reciprocal relation is 
expressed as L rs (d p ) = (-1)KI+KI L sr (d p ) . 
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■ pV 



' 5 AS 

8a v ix) 



(A.14) 



where for a contravariant vector we define Vj_ = (—u^Uu) v v and vft = ft/^ (and similarly 
for covariant vectors). Covariance and positivity further impose the condition that and 
Ljf" should be expressed as If" = Lj_ u^u u (Lj_ > 0) and = L\\ h^ u (Ly > 0) , respectively. 

If we further know the reversible evolutions of thermodynamic variables, [d r (x)] rev , which 
are not accompanied by entropy productions, then the dynamics of the system can be 
determined as 

5AS(t; E X [L 8 }) 



a r (x) = [d r (x)] rev + [d r (x)] irr = [d r (x)] rev + L T 



Sa s (x) 



(A.15) 



For viscoelastic materials, the relevant thermodynamic variables are the following: 



b A 


b A 


c 1 


Pi 


d p 


(d§/dd p ) 


& pV 





P» 


-vP/T 


9tiv 


VhT^/2T 
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-p/T 






e 














where = eu^ u v + is the quasiconservative energy-momentum tensor with the 
quasiconservative stress tensor. The entropy functional is then written as 



AS(t; E X [L S ]) 

= -\l d D y V^gN^x 



( p{v),(pv) pit*"), \ / 

v J V2^v), ^w.w/ \v^{ds/d P(T} ) > 

+ ^ d ^ m V[ er e?) \d„(d~s/dn). 



+ (tre 6 V^ds/dpt) 
( e<H {-l/T)K 



(l\ i\ fc\ ( tre \ 

h h k 

h h h) \v»(ds/dp,)) 



2 (pa) 

Wj 1 „ (pa) f 

2(pv), l 3(pv),(pa), 



S(pa) 

'-l/T)K {pa)) 
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(e, M (-/*/T) 



tre (-1/T) tiK 




(A.16) 



where the coefficient matrices are symmetric and positive semidefinite, and their indices 
are all orthogonal to u M o Note that for this parametrization, the contributions from the 
rotation V^m^j are discarded. Since the matrices must be invariant tensors, we can assume 
that they take the following formic 




9 {liv),{paY 
'2 

9 (fj,u),{pa) 
'3 



where 



ft ft 
*-2 3 



1 *"2 



and 





N h% h v J, W' p h v '° 



(A.17) 
(A.18) 

(A.19) 



are positive semidefinite. The irreversible evolutions 



of thermodynamic variables then become 

1 




5AS 




(5AS/5(tTe) 
5AS/56 



5AS 
5AS 



(A.20) 
(A.21) 
(A.22) 

(A.23) 

(A.24) 
(A.25) 



13 1 



'We actually need to impose the latter condition in the Landau-Lifshitz frame. One can show that if this 
condition is relaxed, the energy-momentum tensor comes to have terms related to a heat flux, which should 
not appear in the Landau-Lifshitz frame. 



14 Note that lW>W e , 



pa) 



2l\e^ (k = 1,2,3). 
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on 



(A.26) 



We now make the following irreducible decompositions of the phenomenological constants 
under the group 0(D) in a local inertial frame: 



■ tr e tr e 



L tree = L 



2 ■ 



r£(i£v — rv L 
Liu = Li ll^y , 

L 90 = L s 3 , 



L nn = M, 



Lj_ M = L± UpU v , L|| M = L|| h^y , 
Then the irreversible evolutions of thermodynamic variables can be written as [13 

iV- 1 [(tr £ )-] ii 



-2L t £\e { ^ ) + (2L t £ t 2 /T)K { ^ ) , 



(A.27) 
(A.28) 
(A.29) 

(A.30) 
(A.31) 



- (L\ £\ + L\ £%) tr e - (L\ i\ + L\ t 3 ) 6 + (L\ £\ + L s 2 t 5 )-trK } (A.32) 



T 
1 



N- 1 [9 ] . rr = - (L s 2 + LI t 2 ) tr e - (LI £1 + L° *|) + (L s 2 + L s 3 t b ) — tr A", (A.33) 



[A], 
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- c ± A± (-u„«a) V„ [ 2 1\ - (2/T) 

+ (£\ tr £ + £ s 5 9 - (1/T) t 6 tr A) /^ A ] , (A.34) 
-qL\\ h uX V^ [2t 2 e^ - (2/T) l\K^ 

+ tr e + 4 - (1/T) t 6 tr A) /^ A ] , (A.35) 

- Vh(-d 2 s/dn 2 )MV fM [t 2 W v e v + £\ W v d„(-fj,/T)] . (A.36) 



Here, in order to evaluate SAS/5p^ , we have used the decomposition of the matrix d 2 s/ dp^dp u 
(negative-definite for each irreducible component) as 

d 2 s 



-c± u»u v - cm hr 



(A.37) 



with positive quantities c± and c\\ . 

If we assume that Lp = c±L± = c\\L\\ and L^ = \/~h(—d 2 s/dh 2 )M are constant, then 
Eqs. ( IA.34j) - flA.36j) are rewritten as 



[PA 



'V"T (d) 



71 



(A.38) 



where the dissipation currents are given by 

rjS = Lp [2£\e {tlv) - (2/T) t\ A w + tr £ + £|0 - (1/T)£ S 6 tr A) V] , (A.39) 
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V 



(d) 



Lnllls^Zlh-dvi-n/T)]. 



(A.40) 



On the other hand, as for the isentropic evolutions, we assume that the evolutions of the 
densities of conserved quantities are given by 



[PuU = -V=9 W« , = -y/=g W« (A.41) 

with the reversible currents of the following form^£ 

t$ = r$ - 2G e^ v) - £ (tr e - a 9) V , (A.42) 

(A.43) 

we set them to be in the most generic 



v 



jJ,U 

(r) 



fW 

-He, 



As for the evolutions of the strains E, 



form: 



N- 1 [e { ^l cv = {2GL"/L i> T)K { , v) , 



'^[(tre)'], 

L J rev 



K'd-iK/T) tiK 



L s 2 L%J \-K'tie + {K,/T)tiK 



(A.44) 
(A.45) 

(A.46) 



Combining Eqs. ( 1A.38I) -( IA.40I) and Eqs. ( 1A.41I) -( )A.43I) . and using the formulas p u 
-gV^u^pv) = y/—g V M (e u^Uy) and n = \J—g V /X (n u^) , we obtain 



^V^eu^,) = = [pj rev + [p v ] i 
v^V^nn,) = ft = [fi] + [ft] 



/=^V(rW+rW), 



(A.47) 
(A.48) 



We thus find that (jO^I - (jOH]l [or Eqs. ( fA~38l) -( fA~4n|l ] and Eqs. (TA~4ll -( fA~43|l can be 
summarized as current conservations: 



V^T' 1 " = , VX = 



where each of the conserved currents, 



(A.49) 



(A.50) 



consists of the convective current {p u u^ = e u^u u or nu^) and the additional current (r M1/ or 
the latter being further decomposed into the reversible and the dissipative currents: 



T - T (r) + T (d) > " - U {v) + ^(d) • 



(A.51) 



15 At the end of this appendix, we comment on how these reversible parts are determined in the entropic 
formulation. 



39 



Furthermore, one can easily show that the evolution equations on E^ v 
[Eqs. flA.30p - flA.33p ] together with the explicit form of the reversible and the dissipative 



currents [Eqs. flA.391) . flA.40p . f)A.42|) . and f)A.43|) ] can be rewritten into the following set of 
equations: 



-(2Ai/T) £ u e 



r (9) 



-(A 2 /T) £ u e, 



( 



■(i/r)7 



■£«(tre)y 
£ u 9 J 
\(1/D) (trr-trr (q) ) / 



2 (£ + 77 



f^ + cr 



(/c + C) 




(A.52) 
(A.53) 

(A.54) 



where 



Ai = 



2D 



A, 



7 




77 



= L ( 





TS TS 



(A.55) 



C = L p 





(A.56) 



(A.57) 



-2 *"4 

;7s «s ps 
-2 c 3 c 5 

( /C' /c \ 
7C= -/C' -/Co 
y-JC JCa J 

Equations flAl^ -f lA~57l) totally agree with Eqs. f l234l) -f l238l) . from which Eqs. f l2T39l -f l2T45l) 
follow, as we see in Sec. 12.31 This is what we promised to show at the beginning of this 
Appendix. 

We close this appendix with a comment on how the reversible evolutions are determined. 
They are actually determined by the requirement that the reversible evolutions do not 
produce entropy and the final form of the total evolutions (reversible ones plus irreversible 
ones) should be given as in Eqs. (1A.52p - flA.57p . As an example, let us consider the irreversible 

(a) 

evolution of e^ u ) and the quantity t^ v ) — tY~L : 



N" 1 [s { ^] iTr = -2Di\e 



1 



2D i\ 



K 



T, 



(p.u) 



r (q) 



— 2£ e^ u ) + 2Lp £ 2 E{ju,) 



{fiu) 5 

2T,~ P t 

P 3 TS~ 

— =: — A, 



(A.58) 
(A.59) 
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By multiplying the first equation by a factor —Lp/I} , the equations can be rewritten with 
a symmetric matrix as 




■(l/r)^>. 



(A.60) 



The second term with a symmetric positive-semidefinite matrix represents irreversible pro- 
cesses with entropy production. Thus, in order for the first term not to produce entropy, we 
need to introduce the reversible part in e^ u ) such that the resulting form can be written with 
an antisymmetric matrix. This consideration determines the reversible evolution uniquely 

as 








(A.61) 



Noting that N 1 iu v \ = £> u Eu v \ , we see that the total evolution is actually given as in 
(IA.52|) . The remaining equations can be obtained in a similar way. 



B. Constitutive equations in rheological models 

The theory of elasticity is based on Hooke's law which states that that stresses are propor- 
tional to strains in elastic materials. On the other hand, the theory of viscous fluids is based 
on Newton's law which states that viscous stresses are proportional to velocity gradients in 
fluids, and is described by the Navier-Stokes equations. However, for more general materials 
these theories are not applicable, and a class of such materials is called viscoelastic materials 
and studied in the area of rheology. The relation between stresses and strains for a given 
material is called the constitutive equations, which play a fundamental role in the study of 
rheology. In this Appendix, we list a few well-known materials with their constitutive equa- 
tions and compare them with the viscoelastic materials discussed in the bulk of the present 
paper. 

Hookean materials 

The simplest constitutive equations constitute Hooke's law. We first assume that, on each 
timeslice H t , every material particle knows its own natural shape described by the reference 
metric hn V , which measures distances in a material when it is free of elastic strains. This 
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metric has the same meaning as the intrinsic metric in the main text, though it is not 
dynamical here (£ u h^ u = 0) . When we discuss nonrelativistic dynamics, we will set it to 
be hfiu = diag(0, 1, . . . , 1) in a laboratory frame, as is taken in standard textbooks (e.g., 
[5]). Although we consider the strain tensor E^ v in the main text, we here assume that 
elastic strains are purely spatial, and only consider the elastic strain tensor defined by 
= (1/2) (hp, - h^u) ■ 

Hooke's law can then be expressed as 

T W = _l C (^)(pa) £ ^ ^u){pa) > q) ^ (B.l) 

where /C^^ is a constant tensor, and ( ) denotes the symmetrization of indices with the 
normalization (()) = ()• For isotropic elastic materials which locally has no preferred 
direction, the coefficient /C^^^ P<J ) can be expressed as the sum of the irreducible components 
hT h pa and (1 /2) (hf* h v ° + h^ h up - (2/D) hT h pa ) , and we have 

Tyt, = -K (tr e) hyt, - 1Q e {fll/) , (B.2) 

where JC and Q are the bulk and the shear modulus, respectively. Relativistic motions of 
such elastic materials in gravitational fields are discussed in, e.g., [T5] . 

Since we are considering the linear approximation in , this stress tensor can also be 
written as 

V, = -K (tire) V - 2Q e-^ , (B.3) 
where tr and are defined by 

tre = h pv , e-^ = 8^- — (tre) . (B.4) 

Then, if we take the nonrelativistic approximation with h^ v = diag(0, 1, . . . , 1) , we reproduce 
the standard Hookean stress tensor 

Hi = -JC {5 kl e M ) 8a - 1Q (e VJ - 1 (5 kl e w )5y) . (B.5) 

The constitutive equations for a Hookean material are schematically represented by a 
spring, as depicted in Fig. HI To understand the diagram, we consider a Hookean material 
in D — 1 spatial dimension. The material can be obtained by connecting in series tiny 
springs with a weight of mass m at each end (see Fig. [5]). Since the actual length between 
two adjacent weights at x = x n = nl (n G Z) and x = x n+ i is given by y/hu(x n ) I , 
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bulk : /C shear : Q 

Figure 4: The bulk part (left) and shear part (right) for a Hookean material, 
a material 

particle > m k 





-10 I 21 

Figure 5: Weights of mass m are connected to the spring with spring constant k . 

and since the natural length is given by \fh\i(x n ) I , the stretch of the spring is given by 
(V^ii — a/ hu) I ~ (su/ \/ h n ) (x n ) , where e n = (1/2) (h n — h n ). Then the equation of 
motion for the weight at x = x n can be written as 

m ai(x n ) = —k I £u (x n „i) + k I £u (x n ) , (B-6) 

where ai(x n ) is the acceleration of the weight at x = x n in the x 1 -direction. Then in the 
continuum limit / — > with e = m/yhxyl and fC = kl kept fixed at finite values, the 
equation becomes 

e ai (x) = [-K (tre{x)) 5\] , (B.7) 

so long as we take a coordinate system in which the intrinsic metric hu is spatially constant. 
If we define the energy density e(x) = m/y/hu I and neglect the difference (e(x) — eo) ai(x) , 
which is of higher orders in e n , we obtain the Euler equation 

e(x)a 1 {x) = -<9it\ (B.8) 

with the stress tensor r\ = —K, (tre) 8\ . This stress tensor coincides with flB.3j) in D — 1 
dimension, and in this sense the left diagram in Fig. H] represents (the bulk part of) the 
constitutive equations for a Hookean material. On the other hand, the right diagram in Fig. 
H] is simply a schematic generalization for the shear part and does not have any physical 
meaning other than the information that the shear part of the stress tensor is given by 
—2Qe^ u ) . 
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Navier- Stokes (Newtonian) fluids 



Newton's viscosity law says that the viscous stress tensor is proportional to velocity gradi- 
ents, and in our notations this can be written as 

T ^ = -C^ v){pu) K p(7 (£(H0») > o) , (B.9) 

because the extrinsic curvature = (1/2) £ u h fJLV can also be expressed as velocity gradi- 
ents, K^u = (1/2) hP hj 7 (V p u a + V a u p ) . In particular, for simple fluids (that do not have 
any specific directions locally) we have 

T^j = -( (tr K) W v - 2r] K^ v) , (B. 10) 

where ( (> 0) and rj (> 0) are the bulk and the shear viscosity, respectively These consti- 
tutive equations can be interpreted as representing the resistance due to the time derivative 
of the induced metric, = (1/2) £ u h IMll , and are schematically represented by a dashpot 
as in Fig. [61 

bulk: o ° shear : O Jrj] O 

Figure 6: The bulk part (left) and the shear part (right) for a Navier-Stokes 
(or Newtonian) fluid. A dashpot yields a viscous stress proportional to the time 
derivative of the induced metric. 

For simple fluids, the reversible part of the stress tensor, rj% , should be proportional to 
W by definition, and we write it as rj"? = P . Then the total stress tensor for simple 
viscous fluids is given by 

T w = T ^ + = Ph' a/ -C (tr K) W v - 2r) K^ v) . (B. 1 1 ) 

Materials with the constitutive equations of this form are called Navier-Stokes (or Newto- 
nian) fluids. 

Kelvin- Voigt materials 

If an elastic material (so that £ u hp U = 0) further obeys Newton's viscosity law, the stress 
tensor is given in the following form: 

V = - (£ tr e + C tr K) V - 2£ £< H - 2rj K {fiu) . (B. 12) 
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Such materials are called Kelvin- Voigt materials and are sometimes used to explain creep 
phenomena in viscoelastic materials. Relativistic motions of such materials are discussed 
in, e.g., [22]. Since Kelvin- Voigt materials have fixed intrinsic metric (0 = £ u h liv = K t 



/IV 



£u£fiu), we have K^ v = £ u e iiv and the stress tensor can be rewritten in the following form: 



Tfn, = - (JC tr e + C £ u tr e) h^ u - 2Q e^ v ) - 2rj £ u e^ u) 



(B.13) 



The constitutive equations for a Kelvin- Voigt material thus can be represented by the di- 
agrams in Fig. [3 Since a spring and a dashpot are connected in parallel in each diagram, 
the total stress is given as the sum of the stress of each component. 



bulk: 



jc 

HaBHSH 



shear : 



-33- 



Figure 7: The bulk part (left) and shear part (right) for a Kelvin- Voigt material. 



Unlike Hookean materials, the stress-strain relation is process-dependent. However, for 
Kelvin- Voigt materials the stress tensor at each moment can be determined only by mea- 
suring the induced metric h^ v and its temporal derivative at the moment, and we do 
not need to know the preceding history of the strains. 

For more generic materials, the stress tensor indeed depends on the whole preceding 
history of the strains. The simplest among such materials are Maxwell materials, described 
below. 

Maxwell materials 

The constitutive equations for a Maxwell material are depicted in Fig. |HJ Since a spring 
bulk : JC shear : Q 

Figure 8: The bulk part (left) and the shear part (right) for a Maxwell material. 

and a dashpot are connected in series, the stress of the spring and the stress of the dashpot 
should be equal. As is already explained, the stress of the spring is given by 

7> = -K (tr e) V - 2Q e {fiu) . (B.14) 
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Recall that the induced metric h^ u measures the actual shape of each material particle (i.e., 
the total length of the diagram), while the intrinsic metric measures the natural shape of 
each material particle (i.e., the length of the dashpot plus the natural length of the spring). 
Thus, the stress of the dashpot, which is proportional to the temporal derivative of (i.e., 
the temporal derivative of the length of the dashpot), is given by 

V = "C (tr K) V - 2r? K (H . (B. 15) 

Since these stresses are equal, from Eqs. (IB. 14|) and (IB. 15|) . we obtain the equations 

tiK= (JC/O tie, K (fM1/) = (G/rj) £ (M ■ (B.16) 

Since K^ v is the temporal derivative of h^ v , K^ u = (1/2) £ u h tll , , these equations describe 
the dynamics of and are called the rheology equations in the main text. Note that the 
structure of a Maxwell material is critically different from that of a Kelvin- Voigt material 
in that the intrinsic metric h^ v of the former is dynamical. 

We should also emphasize that even if we measure the shape of a viscoelastic material, 
, and its derivative K^ v at a given moment, we cannot readily determine the value of the 
stress tensor t^ v at the moment because there is no way to know the values of strains e^ u 
when h^y is dynamical. However, if we observe the evolution of h^ u during a finite interval 
of time, the initial value of can be obtained, and by solving the rheology equations we 
can determine the value of the intrinsic metric h av at each moment. 



Zener materials 



We next consider Zener materials or the standard linear solid model whose constitutive 
equations are given by the diagrams in Fig. [91 This model includes Kelvin- Voigt materials 

shear : 



bulk : JCo 

— mk- 




Figure 9: The bulk part (left) and the shear part (right) for a Zener material. 



and Maxwell materials as limiting cases {K2 = Q2 = and JC± = Q\ = , respectively). 
However, as is clear from Fig. [9j if a Zener material is left intact after an initial deformation, 
it will get back to its original natural shape. In other words, this kind of material does not 
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posses permanent strains unlike Maxwell materials, and in this sense Zener materials are 
said to be solid-like. 

If we want to describe the relativistic dynamics of a Zener material using our theory of 
viscoelasticity, we need to extend the framework, introducing another nondynamical intrinsic 
metric h^l in addition to the original dynamical intrinsic metric h^ v . Here h^o measures 
the natural length of the lower spring in Fig. [91 while measures the length of the dashpot 
plus the natural length of the upper spring. 

If we consider more generic materials, we accordingly should introduce more additional 
intrinsic metrics (dynamical or nondynamical). Such generalizations correspond to consid- 
ering multielement models (such as the generalized Maxwell model) known in the study of 
rheology. In this paper we only consider the cases with a single intrinsic metric, and such 
generalizations will be discussed elsewhere. 



Viscoelastic materials considered in this paper 

As for the rheological model discussed in this paper, we here consider for brevity the case 
when the effects of thermal expansion can be neglected (£2 = Cs = a = I2 = fc' = 0). Then 
the stress tensor and the rheology equations are given by 

V = rj$ -2{Q- 7/2) e^) - (2r/ 3 /T) K {H 

- [(K-U)tTe+{C 6 /T)tTK] V, (B.17) 

P e c i_ g + r ?2 „ fu la , 

J^ u £{fiu) — : — S{jaA H : K {uv) 5 (iJ.loJ 

Ai M 

£ u (tie) = tre+^t^i t rK. (B.19) 

7i 7i 

These equations can be summarized as in Fig. [TO] In particular, one can show that Maxwell's 
original definition is realized if we set JC = 71 — C4 an d G = Ai — r/ 2 (see Sec. 13.3)) . The 
corresponding diagrams are given in Fig. [TTJ The Maxwell model can be obtained if we 
additionally set £e = 773 — , which is the case where the simplified Israel-Stewart model is 
obtained, as shown in Sec. HJ 
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71 (/C-C4) 

bulk : /c+C4 71 (K-a) 



shear : g+m Ai (g-t? 2 ) 



71 (C4-/C)(/C + C4-7l) 



Al (?72-S) (g + 772-Al) 

(S+»72)»;iT 



Cs , (/C-C4)(/C+C4-7i) 
T ' Ci T 



273 _i_ (g-m) (g+>?2-Al) 
T »JiT 



Figure 10: Schematic structure of the bulk part (left) and the shear part (right). 
Note that the contribution of t$ is omitted for simplicity. 



bulk: 



71 (2JC— 71) 

2/C - 7i p 



shear 



Ai (2g-Ai) 

2Q — X\ m T 



Figure 11: A three-element model where a dashpot is connected in parallel with 
a Maxwell material. 
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C. Euler and Gibbs-Duhem relations 

In this appendix we consider the case = P h^ u . Then the variation equation of entropy, 
Eq. (12.24B . is given by 

Ss= -^Sp v - j;5h + ^ SVh + ^ u>V Sg^ 

~2\x Sew ~ ^ (7i tr e + 72 9) *(tr e) 

- ^ A 2 e" Se, - (73 9 + 72 tr 5) SO . (C.l) 
Here, if we consider the variation 5 = w^V^ , we obtain 

y/hu 1 * d^s = u^id^e — ji dn) 

- ^ 2X 1 e<""> u"V p e (H - ^ ( Tl tr e + 72 9) u^(tr e) 

- ^ A 2 ^ w'V^ - ^ (73 + 72 tr e) u^O . (C.2) 
On the other hand, if we consider the variation S = £ u , we obtain 

- ^ 2Ai - ^ (71 tr e + 72 0) u^ M (tre) 

- ^ A 2 e» £ u e, -^( l3 9 + l2 tr 5) u^O . (C.3) 

Subtracting Eq. (10.2)) from Eq. (10.3)) . we obtain the following equation: 

~ , v e- fis + VhP 

s tr K = — trA 

T 

- ^2A X (tr(e 2 K) - (4/D) (tre) tr(eJf) + {A/D 2 ) (tr AT) (tre) 2 ) 

~ "t" 

We can neglect the terms in the second and third lines because they are of higher orders, 
and thus we obtain the equation 

(e + P -Ts- im)XjK = 0. (C.5) 

Since this should hold for any processes in our linear approximations, the following relation 
must hold: 

e + P-Ts- fin = 0. (C.6) 
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X 2 e»e v K^. (C.4) 



This has the same form with the standard Euler relation although the energy density e and 
the entropy density s here contain contributions from the strain tensor = (e^ , , 9) . 
From this and Eq. fIC.lj) . we can derive the Gibbs-Duhem-like equation: 



s 5T + n5fi — 5P 

= -2Ax 6e {H - ( 7l tre + 72 B) 5{tie) - A 2 fe M - (73 9 + 72 tre) 59 . (C.7) 

In the limit where the strains relax completely [E^v — > 0), this reduces to the standard 
Gibbs-Duhem equation for simple fluids. 
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